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Abstract 



A coupled continuum-discrete numerical model is presented to study the synthesis 
of TiC nanosized aggregates during a self-propagating combustion synthesis (SHS) 
process. The overall model describes the transient of the basic mechanisms gov- 
erning the SHS process in a two-dimensional micrometer size geometry system. At 
each time step, the continuum (micrometer scale) model computes the current tem- 
perature field according to the prescribed boundary conditions. The overall system 
domain is discretized with a desired number of uniform computational cells. Each 
cell contains a convenient number of computation particles which represent the ac- 
tual particles mixture. The particle-in-cell (discrete) model maps the temperature 
field from the (continuum) cells to the respective internal particles. Depending on 
the temperature reached by the cell, the titanium particles may undergo a solid- 
liquid transformation. If the distance between the carbon particle and the liquid 
titanium particles is within a certain tolerance they will react and a TiC particle 
will be formed in the cell. Accordingly, the molecular dynamic method will update 
the location of all particles in the cell and the amount of transformation heat ac- 
counted by the cell will be entered into the source term of the (continuum) heat 
conduction equation. The new temperature distribution will progress depending on 
the cells which will time-by-time undergo the chemical reaction. As a demonstration 
of the effectiveness of the overall model some paradigmatic examples are shown. 
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1 Introduction 



The self-propagating high-temperature synthesis process (SHS) [12] is a promis- 
ing method employed for the synthesis of many advanced materials, such as 
ceramic, intermetallics, composites, etc. [13]. Usually the process applies to 
powders mixture which arc conveniently homogenized and pressed in order to 
form a loosely compacted pellet. The SHS exploits the ability of certain mate- 
rials mixture in producing high exothermic and self-sustaining reactions once 
ignited locally or uniformly. Ignition can be made by a laser beam, induction, 
resistance, radiant or a spark source. The exothermic reaction makes the tem- 
perature increase rapidly in the pellet, reaching and surpassing the combustion 
temperature. The resulting product are usually very pure and rather porous 
(about 50% of the theoretical density [13]). The final product composition and 
its morphology depend on [14]: 

i) initial particles size and distribution, shape and purity; 

ii) initial density of the reacting mixture; 

iii) initial temperature of the reacting system; 

iv) size of the sample and reactor configuration; 

v) dilution of the reacting mixture with the final product. 

Thus, the actual combustion process involves simultaneously critical factors 
at both the particle (i.e. i) and the sample scale (i.e. ii-v). 
Typical advantages of the process are: 

• high purity products; 

• low equipment, operation and processing costs; 

• extremely short processing times; 

• possibility of forming unique metastable phases with improved properties as 
a result of the inherent strong non-equilibrium conditions (i.e. steep thermal 
gradients and high heating/coohng rates). 

However, beside the above advantages, the rapidity of SHS and the complex- 
ity of involved concomitant physical and chemical phenomena, makes rather 
difficult the optimization of the process as well as that of the chemistry and 
the morphology of the final products. In addition, the strong exothermicity 
of the reaction may generate a combustion wave which passes through all the 
pellet by igniting abnormally the reaction thus stopping or making unstable 
the reaction. 

It is therefore desirable the development of fiexible predictive simulation tools 

in order to give a significant step forward to the progress of the SHS process. 
In this study we focus our attention to a specific Ti and C particles mixture 
system, in order to obtain TiC particles. According to standard classification, 
this powder system involves a solid-solid reaction in the sense that no gas 
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reactant contributes to the combustion. 

Many studies have been carried out on the SHS synthesis of TiC as model 
material under both microgravity [17,19] and normal gravity [10,18] condi- 
tions, either experimentally [17,19] and theoretically [10,18]. This system is 
particularly convenient because Ti and C powders arc not very expensive. In 
addition, the formation heat of TiC is very large (185 kJ/mol) and the melting 
point of the product is very high (3423 K) which make the reactions initiation 
rather easy and the product formation straightforward. 

The study of this reacting system is complex because of the many processes 
involved e.g. the propagation of the combustion wave in the given granular 
system, the phase changes of the metallic component, the occurrence of the 
chemical reaction, the physical and mechanical interaction among the pow- 
ders and the heat and mass transport in the system which implies strongly 
inhomogeneous properties. 

We propose a discrete approach to model the kinetics of the reaction consid- 
ering the thermal and the chemical evolution of a small portion of the system, 
thereby simulating the actual non uniformity of the mixture properties. 

In the section II, we describe the physical phenomena we take into consid- 
eration in the modelling. Section III describes the numerical methods used 
to reproduce these physics processes at macroscopic and microscopic level 
and their relative coupling. Section IV presents results on the formation of 
nanoparticles during a typical combustion process. 



2 Physical processes and mathematical model 

During the SHS operation, the pellets are rapidly heated. Usually, a ther- 
mal wave precedes the combustion front, thus preheating the powder mixture 
ahead. Usually, the metallic constituent has to melt before reaction occurs 
[15]. Under this condition, the reactions is mainly controlled by the carbon 
diffusion in the liquid Ti. 

The reaction scheme is 

Ti(liq) + C(sol) TiC{sol) (1) 

All the atoms and molecules, and in particular the composites of titanium 
carbide, are subjects to an interaction potential with all the other species and 
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Fig. 1. Enthalpy per unit of mass versus temperature curve for titanium 

at the same time they are thermically pcrtm4:)ed by the Brownian motion. In 
consequence of these competing effects nanoparticles of TiC can be formed. 

In the developed model the following steps are considered: 

a) the solid-liquid transition of titanium; 

b) the chemical reaction between the melted titanium and the graphitic car- 
bon granules; 

c) the interaction between clusters of different kinds of particles (Ti, C and 



d) the formation of nano-sized aggregates of TiC. 

Below we analyze in detail the three fundamental steps: phase change, chemical 
reaction, nanoparticle formation. 

2.1 Heat transfer and phase change 

In the solid-liquid transition the enthalpy of the system has to be taken into 
consideration. In figure 1 the curve of the enthalpy per unit of mass (h) as a 
function of the temperature is shown for titanium. 

The function h{T) includes three continuous parts divided by two discontinuity 
points (i.e. A and B) at the a—f3 and /3-liquid transition, which occur at Ta,/3 — 
1155 K and TM,Ti — 1939 K. During melting or a — /? transition it is necessary 
to supply an extra amount of enthalpy to break the interatomic bonds between 
the atoms of the solid metal (L^ = 295.555 kJ/kg and La,p = 59.247 kJ/kg 
respectively) . 

The enthalpy versus temperature function for a mass m is defined as 



TiC); 



n{T) = mh{T). 



(2) 
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The heat conduction equation for the system is given by [9] : 



pC^^VkVT + S, (3) 

where T is temperature, S the energy source, p the density, C the specific heat 
at constant pressure and k the thermal conductivity. 

However, equation (3) it is not convenient in this form when 7i{T) is dis- 
continuous. In this condition the presence of a heat fiux does not modify the 
temperature in the volume V where phase change is taking place and the eq. 
(3) can be replaced more conveniently by: 



if hp^f q (4) 
at Jv Jdv 

where q is the specific heat fiux (J/m^s) across its surface given by the standard 
Fourier law 



q = kVT. (5) 



2.2 Chemical reaction 



As experiments suggest we assume that the chemical reaction occurs when 
the Ti surrounding the C particles has melted, while the carbon remains solid 
since its melting point, Tm,c = 5000 K, is far higher temperature. Similarly 
the TiC being created is solid since its melting temperature is Tm^tic = 3423 
K. 

Due to the exothermicity of the reaction, the heat source in eq. (3) is evaluated 
as 



•S* = niTic^Hr, (6) 

where niTic is the mass of the TiC product and AHr indicates the heat gen- 
erated per unit of mass of product(Ai7j. = 3.08 • 10^ kJ/kg [7]). 
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2.3 Nanoclusters interactions 



After the reaction has taken place, the temperature of the system is locally 
increased and the TiC product becomes dispersed in the liquid Ti phase. Since 
the atomistic simulation of Ti, C and TiC is computationally prohibitively 
expensive, we consider the dynamics of clusters made of a large number of 
atoms. The behavior of the clusters attemps to reproduce on the average that 
of the single particles relatively to the phenomena we study. 
The motion of the cluster is partly deterministic (due to the interaction with 
the other clusters) and partly stochastic, because of their Brownian motion. 
For the particles / and m the potential energy has the form [3] 



a 
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where r^^ = r; — and r — |r|. The parameter e defines the strength of the 
pair interaction and a defines a length scale [3]. The interaction is repulsive 
at close distance, and attractive at larger distance. A cut off is assumed, for 
numerical convenience, beyond a limiting separation Vc- The corresponding 
Lennard- Jones force is F = — V^(r) and has the expression 
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The detailed motion of the clusters is described by the Langevin equation [8] 

m^ = -av + F*(t) (9) 

where m is the cluster mass, v the velocity, a a parameter depending on the 
liquid [8] and F*{t) the Brownian random force acting on the cluster. 

The motion of these clusters, and in particular that of the TiC clusters, can 
generate a phenomenon of coalescence in the form of very small aggregates of 
nano-metric sizes. Thus, starting from micro-sized reacting particles of carbon, 
we obtain nano-sized aggregates of TiC. 



3 Numerical methods 



The simulation of the processes described above and the discretization of the 
various mathematical models used to describe them requires to handle multiple 



6 



Discrete 



Properties of 
particles 



9 



Nanoscale 




Heat Transfer 
+ 

Reactions 



Temperature 
of ceiis 



Position and veiocity 
of particies 



PIC 

{X,V}iP 

I t 



{k, C, ?p 



Properties of 
ceiis 



Mesoscale 



'{X,V}2P 



r 



Temperature 
of particles 



^2" 



MD 



Fig. 2. Coupling of the continuous and the particle models 

length and time scales. The Brownian motion and the chemical reactions are 
much slower than the evolution of the temperature. The processes at the 
particle level (chemical reactions, motion) require a spatial resolution at the 
nanometer range, while the temperature field is characterized by micrometer 
scales and the overall system is on the centimeter scale. 



To handle these multiple scales we have decided to focus our attention at 
the mesoscopic level, describing the evolution of only a portion of the com- 
plete system. Our domain is on the micrometer range and is treated at that 
level using a continuum model discretized with a finite volume (FV) method. 
However, we retain the physics at the nanometer range by resorting to two 
additional methods: we describe the nanometer physics using computational 
particles. The interactions among particles at the nanometer range are treated 
with the molecular dynamics (MD) approach, and the coupling between meso- 
scopic (micrometer) scales and nanometer scales is treated with the particle 
in cell (PIC) method. Figure 2 summarizes the approach. 



This approach has the advantage to use the numerical approach most suitable 
for each phenomenon we consider but requires to handle multiple methods 
and algorithms and their mutual interaction. To this end we have applied a 
modular approach based on the object-oriented software paradigm. 
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3.1 Mesoscopic model 



The heat transfer in a system with phase changes is described by equation (3) 
coupled to equation (5). It is a continuum model that involves a local energy 
balance. At each point of the system the enthalpy content, which depends on 
the local temperature and on the state of aggregation of the matter, is known. 
When a chemical reaction event occurs, the amount AHr of energy released is 
added to the local enthalpy content as a source term in eq. (3). Consequently 
the local temperature is increased. The TiC particles generated become a part 
of the system. 

The integration of (3) and (5) is the obtained using a finite volume method on 
a two dimensional cartesian domain Q of sizes l^, ly with a boundary indicated 
by r. This system is discretized using a two dimensional uniform grid in the 
direction x and y with cells of sizes Ax and Ay. We also defined a third 
dimension of size l^ which is used to dimensionalize all quantities as in the 
physical system. The time is discretized in time steps At. 
The discretized forms of eq. (4) and (3) are applied to cells with or without 
phase change respectively. To avoid numerical instabilities we use an implicit 
scheme based on the Euler algorithm. 

• cells with no phase change 



CI Pi ' ' = K Ti — + (10) 



rnn + l rpn rpn+l _j_ rpn + 1 ryrpn+l 

rpn+l _|_ rpn + 1 c^rpn+l 

I un i,j I qn+1 



• cells with phase change 



AH^j = A,Atq^, (11) 



being Aij — 2lz{Ax + Ay) the perimeter of the cell in which phase change is 
taking place. The heat flux is calculated by the discretized form of eq. (5): 



„n ijn* 



rpn rpn rpn rpn 

in* ij I un* ij i+'^J /lo^ 

•iij-l^i-lj +l^i+l,j U^J 

rpn T^TI rpn rpji 

where kf*i j is the thermal conductivity at the interface between the i, j and 
the 1 — 1,3 cells. Its value is given by 

un I un 

t€i, - (13) 
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Analogous expressions are valid for all other interfaces. 

Upon computing the coefficients kij at n the time level , the eqs. (10)-(13) 
become a linear system which is solved using the GMRES algorithm. 

A new algorithm has been implemented to solve the problems arising from 
the non-linearity and discontinuity of the enthalpy function versus tempera- 
ture, 7i{T). In order to illustrate the approach, we consider for convenience 
solidification first, and then melting for a given cell of the system. 



3.1.1 Solidification 

For temperatures above TM,Th we use eq. (10) to advance the temperature T^" 
for the cell ij at the temporal step n , the corresponding enthalpy content for 
a generic cell is given by 



m = n{TJ^A. (14) 



13 ^3 

We indicate with Hi and Hj {Hi > Hi) the enthalpies of a cell of titanium 
corresponding to the beginning and the end of the solidification respectively. 

Solidification starts when for tt. = n we have 

hI>Hl 

< TM,Ti (15) 
Hf+^ < Hi 



and at subsequent time levels m > n -|- 1 , the following condition holds: 

(16) 



> Hi, 



In that case, eq. (10) is no longer appropriate and the enthalpy and the tem- 
perature are advanced using: 

T'm ^ ^ 

ij — J-M,Ti- 
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Finally, the solidification process is completed at the time step m — m, when 
the following condition holds: 

H^+^ <Hi<H^ (18) 



In this last instance, the temperature is advanced as 

T^^^ = n-\Tf^^'). (19) 



Later, soUdification no longer occurs in the cell and, for time steps r > m + 1, 
the temperature is again computed using eq. (10) and the enthalpy is given 
by 

Hi^^nirr). (20) 



3.1.2 Melting 

In the case of melting, the system starts at a temperature below Tjv/.Ti and eq. 
(10) governs the evolution of the temperature. Again the enthalpy is computed 
from 

= n{T^^) (21) 



The process of solidification starts at the time step u = u, when the following 
conditions are verified 

' > TM,T^ (22) 



Subsequently, ion: v > u + 1, enthalpy and temperature are calculated using 

H-j = H^tj + AijAt EL«+i 4 (23) 

ij - J-M,Ti- 



Solidification ends at the time step v — when the condition below is verified 

hI<Hl< Hf;' (24) 
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At the time step v + 1, the temperature is advanced as prescribed by conser- 
vation of energy: 

Tg+^ = 7i-^(7;f^). (25) 

Finally when melting is completed, ior z > v + 1 the temperature is advanced 
using eq. (10) and the enthalpy is given by 

HI, = (26) 



3.2 Coupling of mesoscopic and nanometer scale through PIC 

The mesoscopic subdivision of the system is made using the grid of cells de- 
scribed above. The cells are used to treat the heat transfer in a melting or 
solidifying continuum media. 

To study the nanometer scale physics and the detail of the phenomena at par- 
ticle (i.e. nanometer) level, it is necessary to use a discrete approach. Thus, 
the Particle In Cell method (PIC) [2] has been adopted. This implies that the 
system is subdivided in computational particles which physically represent a 
cluster of atoms. Therefore the granules of carbon correspond to a set of com- 
putational particles with the properties of the carbon atoms. Similarly the Ti 
and the TiC arc represented by computational particles with the correspond- 
ing properties. Then the cell properties are computed based on the particles 
properties which change in time. This approach is particularly convenient in 
the study of powder systems when the mass transfer and the heat transport 
are closely coupled. 

The two levels of description of the system permit us to map not only the ther- 
mophysical properties (C, p and k) from one level to the other, but also to map 
the field variables such as temperature and enthalpy. Thus the temperature 
and the enthalpy of a particle can be computed by knowing the corresponding 
values in the cells and viceversa. 

The approach followed is summarized in Fig. 2. The system is characterized 
by Np particles, each of volume Vp, and concurrently the same system can also 
be represented by cells of volume Vc- 

Assuming an equal number of particles per cell, Npc, the initial volume of each 
particle is given by 



Prom particles to cells 


Prom cells to particles 




C 


" T.,VpW,p 




^ 12p Cp^pVpWcp 

Z^p Vpyvcp 





Table 1 

Properties and physical quantities calculated using PIC 

By introducing a system of logical coordinates ^ i — ^-nd rj — for 
a generic cell c and particle p we define the assignment function W as the 
function resulting from the tensor product between h — spline functions of 
first order in each direction: 



W{ic - 4' nc - lip) = 6i(^c - 4)^1 (^c - (28) 

where (^p, r]p) is the position of the particle, (^cj Vc) the position of the cell 
center and hi is given by 

otherwise. 



This function which will be denoted by Wcp for shorthand, gives the contri- 
bution of a generic property of the particle p to the cell c as a function of the 
relative positions. 

Similarly, it is possible to perform the opposite operation consisting in map- 
ping a property (or a physical quantity) of the cells to a particle. 

The properties and the physical quantities that are calculated with this method 
are summarized in Table I. By knowing the temperature in each cell, at each 
time step it is possible to single out those where the reaction can take place. 
By recalling the basic assumption, e.g. the titanium must be melted before 
the reaction can take place, the condition that must be verified is 



> TM,T^. (30) 

Another necessary condition in order to have a reaction in the cells which fulfil 
the condition (30), is that each cell contains at least one particle of titanium 
and one of carbon. If more than one couple is present, then the reaction 
occurs first between those particles which have the shortest distance. These 
two particles are thus replaced by two new ones. If we indicate with Ti the 
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particle of titanium and with C that of carbon, the number of moles of the 
new TiC particle TiC being created is 



< 



riTic^nTi if nTi<nc 
riTiC = nc if nn > nc 
nric^riTi if nn^nc 



(31) 



whereas its volume is 



VtiC — 



riTiciPMn + PMc) 



(32) 



PTiC 



where PMti and PMc are the molecular weights of titanium and carbon. 
For sake of mass balance a new particle must be introduced into the system. 
This can be either C or Ti particle. Its volume must fulfils the conditions (33) 



If the condition n^i = nc is verified, only a particle of TiC is created. 

The amount of energy developed during the formation of the TiC particle is 

calculated using eq. (6). 

3.3 Molecular Dynamics treatment of nanometer scale particles 

The interaction between difi^erent particles in the system and among the atoms 
that compose them is an aspect which plays a very important role in the 
phenomenon of the formation of nanoparticles. This interaction is reproduced 
by the Lennard- Jones potential. This is only a first approximation of the 
real situation, because it usually holds for atoms or molecules and not for 
clusters. Indeed, our computational particles represent cluster of atoms that, 
on average, we assume interacting with a potential that includes a repulsive 
and an attractive part as that of Lennard- Jones. In calculating the potential 
we consider the distance between the center of the clusters. 
The equations of the motion for a generic computational particle are 




(33) 




(34) 



(35) 
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where x, v and F are position, velocity and force for a particle of mass m. The 
equations of motion are solved using the Verlet algorithm [1]. To simulate the 
stochastic behaviour due to the Brownian motion of the particle we use the 
Andersen thermostat approach [1]. We consider each cell as interacting with 
an heat bath which sets its temperature at the value given by the continuum 
model described above while, at the same time, it permits this subsystem to 
access all the energy shells corresponding to this temperature according to 
their Boltzmann weight [1]. The coupling of the system to the bath is rep- 
resented by stochastic velocities which are occasionally randomly assigned to 
the computational particles, using a collision probability. The strength of this 
coupling depends on the frequency z/ of stochastic collisions between particles. 
The values of the accessible velocities are those of the Maxwell-Boltzmann 
distribution [8] with variance 



a=J^, (36) 
V m 

where ks is the Boltzmann constant and Tp is the temperature of the particle. 



4 Results of the simulations 



4-1 Validation test - ID analytical benchmark 



At first, we performed a comparison between the results of the present model 
with an analytical solution for a system undergoing a phase change [9]. The 
simulation has been carried out on a one-dimensional domain of length I — 
0.1 m made of titanium for a time tf — 150 s with the following initial and 
boundary conditions 



< To" = 2500 Vn (37) 



In fig. 3 we show the time evolution of the liquid-solid interface. As it is 
possibile to see the agreement with the analytical solution is good. We have 
conducted a convergence study to verify the correctness of the implementation 
of the heat transport algorithm described above. 
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Fig. 4. Computed enthalpy for the central cell of a solidifying system made of 
titanium 



4-2 Validation test - 2D phase change benchmark 



We carried out a simulation for a two-dimensional system made of titanium 
with sizes Ix = ly — 0.015 m for a time tf — 5 s. The initial and boundary 
conditions were 



300 + 3000sm (^x^) Vi 



T^ = 3QQ Vn. 



(38) 



The enthalpy versus the temperature for the central cell is shown in fig. 4. The 
result proves that the phase change happens correctly and with the correct 
enthalpy change corresponding to the latent heat. 
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4-3 Formation kinetics of TiC nanoparticles 



In this case the system considered is composed by titanium surrounding three 
granules of graphite, two with a diameter 1.5 iim and the other with diameter 
1 /im. The size of each side of the domain is Ix — ly — 10 //m. 
The properties of the materials, and in particular the specific heat and the 
thermal conductivity, depend on the temperature and are obtained from [4,5]. 
The density is assumed as constant for both materials and equal to that at 
300 K. This means that the effects of the volume variation of the species 
during the process are neglected. 

A generic heat source is introduced to simulate the melting of titanium. The 
thermal behavior of the system assumes the following initial and boundary 
conditions 



The domain is subdivided by using A^^, and Ny cells along x and y directions 
respectively and with Np^ particles per cell. The simulation time is t/, the 
numbers of time steps in the heat equation and in the motion equations are 
nth and respectively. Subcycling of the heat equation is used to handle 
its much faster scale, compared with the scale of the particle motion. For the 
simulation shown below, the following values were used: = 10, Ny = 10, 
TVpc = 9, = 1 • 10-^ s, nth = 1500 and = 50000. The initial velocity of 
all particles is zero, allowing the Andersen thermostat to establish the proper 
thermal equilibrium. The collision frequency used in the Andersen thermostat 
is u — 2b kHz. 

The time evolution of the reacting system is shown in fig. 5. The evolution of 
the combustion process is evident: TiC shells of particles are first formed at the 
interface between graphite and liquid Ti, as demonstrated in the experiments. 
After this stage, the thermal agitation of the particles breaks the shells and 
permits that the reaction develops further. The final aggregates of particles 
are the end product of the reaction: TiC nanoparticles. 

Figure 6 shows the evolution of the concentration of the various chemical 
species for the same simulation. 




(39) 
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Fig. 5. Evolution of the system 




Fig. 6. Evolution of the concentration of the chemical species 
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5 Conclusions 



A new coupled mesoscopic-nanometer scale model has been developed to sim- 
ulate the SHS process. The approach combines the solution of PDEs for the 

continuum media on the mesoscopic level with the discrete PIC-MD tech- 
niques for the nanometer scale. This model permits the description of thermal 
phenomena on the mesoscale and physical and chemical interaction on the 
nanoscale of both reagents and products species. 

The overall model offers a unique opportunity to follow the evolution of the 

nanostructure during the SHS process. 

An interesting output of the meso/nano model is the concentration history of 
chemical species which are relatively easy accessible by actual experiments. 
The model has been validated against an analytical solution. The PIC-MD 
model requires direct experimentation in order to determine the most suitable 
LJ parameters and the collision frequency value for more realistic application 
of the developed model. 
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